# This is the R code used to do data analyses in the article "Terrorism, Counterterrorism Aid, and Foreign Direct Investment." 

setwd("Replication files_FPA")

data <- read.csv("FPA_terrorism&FDI_data.csv") #read data

## Table 1 ##
m1 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount+factor(country), data=subset(data, year>1986))   #Model 1
m2 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount+lnctaid+factor(country), data=subset(data, year>1986))  #Model 2
m3 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount*lnctaid+factor(country), data=subset(data, year>1986))  #Model 3
m4 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount+lnUSaid+factor(country), data=subset(data, year>1986))    #Model 4
m5 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount*lnUSaid+factor(country), data=subset(data, year>1986)) #Model 5

## Table 2 ##
m2.1 <- lm(I(FDI.GDP*100)~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount*lnctaid+factor(country), data=subset(data, year>1986)) # Robust Model 1
m2.2 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount*lncaid+factor(country), data=subset(data, year>1989))  #Robustness Model 2
m2.3 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lngtd*lnctaid+factor(country), data=subset(data3, year>1989)) #Robustness Model 3
m2.4 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lngtd*lncaid+factor(country), data=subset(data3, year>1989)) #Robustness Model 4
m2.5 <- lm(lnfdi~lnGDP+I(GDPpc/100)+growth+trade+sexchange+polity2+durable+population+I(sum/1000000000)+oil+inwar+atopally+sanction+lncount*lnctaid+factor(country), data=subset(data3, year>1986)) #Robustness Model 5

## Figure 2: Marginal effect of terrorism ##
quantile(data$lnctaid, probs=c(0.001,0.01,0.05,0.1,0.2,0.3,0.35, 0.4, 0.5, 0.8, 0.81, 0.85, 0.9, 0.95, 0.99, 0.995, 0.999, 1), na.rm=TRUE)
x1 <- seq(from=-1, to=17.41, length=1000) # for graph
gamma0 <- coef(m3)["lncount"]
gamma1 <- coef(m3)["lncount:lnctaid"]
marginal.effect <- gamma0 + gamma1*x1

var.gamma0 <- (summary(m3)$coef[,2]["lncount"])^2
var.gamma1 <- (summary(m3)$coef[,2]["lncount:lnctaid"])^2
cov.g0g1 <- (vcov(m3)[15,16])*sqrt(var.gamma0)*sqrt(var.gamma1)
var.mar.eff <- var.gamma0 + var.gamma1*x1^2 + 2*x1*cov.g0g1
se.count <- sqrt(var.mar.eff)

ci <- cbind(marginal.effect-1.959964*se.count, marginal.effect, marginal.effect+1.959964*se.count)

par(mar=c(4,4,0.5,1))
matplot(x1, ci, lty=c(2,1,2), ylim=c(-2, 5), type="l", col="black", lwd=2, main="", xlab="", ylab="", cex.axis=1.2)
abline(h=0, col="gray30")
rug(data$lnctaid)
rug(seq(-0.5,0.5, length=1000))
mtext(side=1, "US Counterterrorism Aid (logged)", line=2, cex=1.2)
mtext(side=2, "Marginal Effect of Terrorism", line=2.2, cex=1.2)
legend("topleft", c("Marginal effect", "95% confidence intervals"), cex=1.2, lty=c(1,2))

